home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / BESY.PAS next >
Pascal/Delphi Source File  |  1985-04-03  |  2KB  |  91 lines
  1. program besy;        { -> 348 }
  2. { evaluation of Bessel function of the second kind }
  3.  
  4. var    x,ordr    : real;
  5.     done    : boolean;
  6.  
  7. external procedure cls;
  8.  
  9. function bessy(x,n: real): real;
  10. { cylindical bessel function of the second kind }
  11. const    small    = 1.0E-8;
  12.     euler    = 0.57721566;
  13.     pi    = 3.1415926;
  14.     pi2    = 0.63661977;    { 2/pi }
  15. var    j    : integer;
  16.  
  17.     x2,sum,sum2,t,t2,
  18.     ts,term,xx,y0,y1,
  19.     ya,yb,yc,ans,a,b,
  20.     sina,cosa        : real;
  21.  
  22. begin            { function bessy }
  23.   if x<12 then
  24.     begin
  25.       xx:=0.5*x;
  26.       x2:=xx*xx;
  27.       t:=ln(xx)+euler;
  28.       sum:=0.0;
  29.       term:=t;
  30.       y0:=t;
  31.       j:=0;
  32.       repeat
  33.     j:=j+1;
  34.     if j<>1 then sum:=sum+1/(j-1);
  35.     ts:=t-sum;
  36.     term:=-x2*term/(j*j)*(1-1/(j*ts));
  37.     y0:=y0+term
  38.       until abs(term)<small;
  39.       term:=xx*(t-0.5);
  40.       sum:=0.0;
  41.       y1:=term;
  42.       j:=1;
  43.       repeat
  44.     j:=j+1;
  45.     sum:=sum+1/(j-1);
  46.     ts:=t-sum;
  47.     term:=(-x2*term)/(j*(j-1))*((ts-0.5/j)/(ts+0.5/(j-1)));
  48.     y1:=y1+term
  49.       until abs(term)<small;
  50.       y0:=pi2*y0;
  51.       y1:=pi2*(y1-1/x);
  52.       if n=0.0 then ans:=y0
  53.       else if n=1.0 then ans:=y1
  54.       else
  55.     begin        { find y by recursion }
  56.       ts:=2.0/x;
  57.       ya:=y0;
  58.       yb:=y1;
  59.       for j:=2 to trunc(n+0.01) do
  60.         begin
  61.           yc:=ts*(j-1)*yb-ya;
  62.           ya:=yb;
  63.           yb:=yc
  64.         end;
  65.       ans:=yc
  66.     end;
  67.       bessy:=ans;
  68.     end        { x<12 }
  69.   else        { x>11, asymtotic expansion }
  70.     bessy:=sqrt(2/(pi*x))*sin(x-pi/4-n*pi/2)
  71. end;    { function bessy }
  72.  
  73. begin
  74.   cls;
  75.   done:=false;
  76.   writeln;
  77.   repeat
  78.     write('Order? ');
  79.     readln(ordr);
  80.     if ordr<0.0 then done:=true
  81.     else
  82.       begin
  83.     repeat
  84.       write('Arg? ');
  85.       readln(x)
  86.     until x>=0.0;
  87.       writeln('Y Bessel is ',bessy(x,ordr))
  88.     end        { if }
  89.   until done
  90. end.
  91.